T94/007 

EOTmat/ 0402043 

Finite-lattice extrapolations for a 
Haldane gap antiferromagnet 



O. Golinelli*, Th. Jolicceur^ and R. Lacaze^ 
Service de Physique Theorique, CEA Saclay 
F-91191 Gif-sur-Yvette, France 



Abstract 



We present results of exact diagonalizations of the isotropic antiferromagnetic 
S = 1 Heisenberg chain by the Lanczos method, for finite rings of up to 
N = 22 sites. The Haldane gap G(N) and the ground state energy per 
site e(N) converge, with increasing N, faster than a power law. By VBS 
and Shanks transformations, the extrapolated values are G(oo) = 0.41049(2) 
and e(oo) = —1.401485(2). The spin-spin correlation function is well fit by 
exp(— rj£)l^fr with £ = 6.2. 

PACS number: 75.10.Jm 



Typeset using REVTeX 



*email: golinel, thierry and lacaze@amoco.saclay.cea.fr 
tCNRS Research Fellow. 

1 



I. INTRODUCTION 



A great variety of magnetic phenomena can be understood by the study of classical 
spin systems. However, we know that there are some surprises from quantum mechanics in 
one dimensional spin systems. Haldane has conjectured |lj that the properties of the one- 
dimensional Heisenberg antiferromagnet are qualitatively different depending on whether 
the spin is integer or half-integer. This intriguing argument applies notably to the simple 
prototypical antiferromagnetic (AF) S = 1 spin chain which, according to Haldane, should 
exhibit a gap (G) towards spin excitations. This conjecture has been checked experimen- 
tally, in particular with the compound NENP, for which inelastic neutron scattering and 
susceptibility measurements have clearly shown the existence of a spin gap @ . 

The Heisenberg AF S = 1 spin chain has been studied in many numerical works. In 
1973, ten years before the Haldane conjecture, De Neef || computed the specific heat by 
exact diagonalizations of the Hamiltonian for chain lengths up to iV = 8. In 1977, Blote 0] 
diagonalized chains up to N = 10. In 1982, Botet and Jullien ||, with diagonalizations up 
to iV = 12, obtained evidence for a non-vanishing gap in the thermodynamic limit. Their 
value for the gap was rather imprecise (G ~ 0.25 J). After this work, other authors used 
exact diagonalizations with chains of increasing length: in 1984, Glaus and Schneider 
and Parkinson and Bonner [7| with N = 14, in 1987, Moreo H with N = 16 and in 1990, 



Lin || with iV = 18, a length that has also been reached by Haas et al. |10| and Takahashi 
jl 1| . This growth is almost linear: two more spins every three years. In fact, the numerical 
complexity of an exact diagonalization is 3^. The exponential growth of computer power is 
not sufficient to explain these results and much progress has been done in the algorithms. 
However it is clear that it is very difficult to continue along this line of study. 



With Monte-Carlo methods |ll|-|15|, longer chains can easily be studied (for example 



N = 64), but the results have statistical as well as systematic errors (the latter being much 
more troublesome). In 1985, Nightingale and Blote [16 obtained a very precise estimate of 
the Haldane gap G = 0.41 by use of a stochastic implementation of direct iteration. In this 
case there is a systematic bias, caused by the finite number of "walkers", that is difficult to 
control. 

Real-space renormalization-group methods have also been applied to spin chains. The 
first works were quite disappointing because of large systematic errors, but in 1988 Lin and 



Pan [|17]], gave the very precise value G = 0.4097(5). In 1992, White and Huse [jn| improved 
this method and found G = 0.41050(2) and a ground state energy with a precision of 10 -12 . 

In this work we have used the Lanczos technique on the longest possible chain we could 
reach which is N = 22 spins and then we have applied powerful extrapolation techniques. 
In this strategy, the only source of error is due to the extrapolation technique while the 
finite-lattice numbers are limited essentially by machine accuracy. 

In Section II, we explain the numerical method (in particular the importance of the 
symmetries and the quantum numbers of the Haldane triplet). The programming techniques 
useful for a chain length N = 22 are described in an Appendix. In Section III, we explain 
our extrapolation method: the Shanks and VBS transformations and how we quote errors. 
In Section IV, we apply our strategy to the Haldane gap and the ground state energy. In 
Section V, we compare our results with those of other authors. In particular, the precision 



for the gap is similar to that of White and Huse [18] with compatible results. In Section 
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VI, our results for the correlations functions are presented. They are well described by a 
correlation length £ = 6.2. Section VII presents our conclusions. 



II. NUMERICAL METHOD 

The Hamiltonian for a chain of length N is: 

N 

H = J2S l -S i+1 , (1) 
1=1 

where the Si are the quantum spin-1 operators. The exchange constant is positive (J = 1) 
in the antiferromagnetic case. The boundary conditions are periodic (N+l = 1) and the 
length iV is even to avoid frustration. 

We have computed by exact diagonalization of the hermitian matrix H, the ground state 
|0), its energy Eq, its spin-spin correlation functions and the energy of the first excitation E\, 
for finite lengths up to iV = 22 spins. The Haldane gap G(N) is the difference between E 
and Ei, the two lowest eigenvalues. The extrapolation method, which gives an estimation 



of G(oo) is explained in the next section. We have used the standard Lanczos method [|19 
which is well adapted for this problem: the matrix is very sparse and only a few extreme 
eigenvalues are needed. 

The size of the matrix H is3 N x 3 N . With the symmetries of hamiltonian, H is block di- 
agonal and only the two interesting blocks, with a size smaller than 3^, must be diagonalized. 
The symmetries of the lattice are the translational invariance T, and the left-right reflection 
Lr (Lr transforms the wave vector k to — k; so it reduces the size of blocks only for k = or 
7r). The spin symmetries are the global rotation, S = J^i^i (it seems difficult to implement 
this symmetry, and in practice, only a sub-group of SU (2) is used). The matrix elements are 
computed in the z-axis basis: {\si, S2, ■ ■ ■ , sn)} where Sj = —1, or 1, are the eigenvalues of 
S-. In this basis, T\si, . . . , sn) = |s2, • • • , $n, si) and Lr\{si}i) = \{sN-i}i)- The spin sym- 
metries that can be implemented easily are: S z , the magnetization along the z-axis (which 
is diagonal in the z-axis basis) and the 7r rotation around the x-axis, R x = exp(iirS x ). In 
this basis, R x \{si}i) = |{ — Sj}i) and the action of R x is a flip of all the spins. So, R x maps 
the subspace S z = m onto S z = —m and reduces the size of blocks only for m = 0. 

We will now explain which blocks contain the ground state |0), and the first excitation 
Because of the SU (2) symmetry, the eigenvectors can be labeled by the quantum numbers 
j and m. Each energy level has a degeneracy 2j + 1 and a representative member in the 
subspace m — 0. On the other hand, the subspace m — 1 contains no singlet j = 0. With 



the help of general arguments ||20|| , the ground state |0) of an antiferromagnetic model is a 
singlet j = 0, but the first excitation can have j = or 1. The full diagonalization of H 
for short chains shows that the first excitation has j = 1: the Haldane triplet. We denote 
it with m = —1,0 or 1. The other quantum numbers are obtained by using the 

transformation [EU 

U = ex V Ln £ sA. (2) 

V j even / 
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which is diagonal in the z-axis basis. The interest of U is that the non-diagonal elements 



of UHU^ 1 are or —1. The Perron-Frobenius theorem [21] can then be applied in each 
subspace S z = m. For m = or 1, it follows that the components of U\0) and U\l, 1) are 
strictly positive (in the z-axis basis). In a subspace S z = m, T ■ U = (—l) m U ■ T. So |0) has 
the wave vector k = and |1, 1) (and so the entire Haldane triplet \l,m)) has k = ir. The 
left-right reflection Lr commutes with U; |0) and |1, 1) (therefore all the triplet) are even 
for Lr. The spin rotation R x commutes with U; |0) is even for R x . But R x \l, 1) = ±|1, — 1) 
and only |1,0) is a eigenvector of R x . In a triplet, the eigenvalues of S x are —1,0 and 1; 
so those of R x are —1, 1 and —1. The eigenvectors |1, 1) ± 1 1 , — X ^ match ±1; necessarily, 
the eigenvalue of |1,0) for R x is —1. To summarize, the ground state |0) is in the subspace 
(S z = 0, k = 0, Lr = +1,R X = +1) and one representative of the Haldane triplet is the 
ground state of the subspace (S z = 0, k = it, Lr = +1, R x = —1). 

The advantage of these symmetries is the reduction of the size of the matrix. The size 
of the subspace S z = m is N\/(n + \ n \ n_!) with n + + n Q + n_ = N and n + — n_ = m. 
When N is large, 

^ = 0)--^=. (3, 

The translation T reduces this size by a factor iV (asymptotically when is large), the 
left-right reflection Lr by 2 (N large) and the 7r-rotation R x by 2 (N large). For N = 22, 
the size is reduced by 851 (1% better than the asymptotic formula) and it is equivalent to 
N = 16 without symmetries. 

Certain methods are well adapted to obtain the ground state of a large, unstructured 
sparse and symmetric (or hermitian) matrix, for example the conjugate gradient |22| and 
the Lanczos |T5[ method. In these iterative methods, by starting with an arbitrarily vector 



Vo, the matrix H acts only in matrix- vector multiplications and remains sparse: the vector 
V n is a linear combination of H ■ V n -i and the previous Vf's (i < n — 1). Then V n is in the 
subspace fC n = span(Vo, H • V , . . . , H n ■ V ). From a theoretical point of view, the Lanczos 
method builds an orthogonal basis of fC n and the projection of H on K n is a tridiagonal 
matrix n x n. After n iterations, the ground state is approximated by the vector of K, n 
which minimizes the Rayleigh quotient R(V) = (V\H\V) / (V\V). So, by construction, it 
is the fastest convergent method among these one using H ■ V multiplications. Because 
computers have a finite precision, the orthogonality of the Vi tends to be lost after many 
iterations and it is difficult to obtain many eigenvalues. However, as we want only the 
ground states of some blocks, we used the Lanczos method. For N = 22, only 55 iterations 
(or matrix-vector multiplications) are needed to obtain eigenvalues with a precision which 
can not be improved by more iterations. Details on the programming techniques are given 
in Appendix. 



III. NUMERICAL RESULTS AND EXTRAPOLATION METHOD 

The numerical values are given in Table | with 12 digits after the decimal point, for 
periodic chains with length up to = 22. We have no direct means to estimate the 
precision. Some direct iterations have been made with the eigenvector obtained by the 
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Lanczos method and the precision is estimated at better than 10~ n for the ground state 
energy E and the first excitation Ex- The sizes N = 2 and 4 are added because we will 
see that they have, surprisingly, a good behavior with respect to the extrapolation method. 
Results up to iV = 18 have been published by other authors (see the caption of Table |). 

The gap G(N) — Ex — Eq and the ground state energy per site e(N) = E /N have still 
not converged. To obtain a good estimate of their thermodynamic limits, the convergence 
must be accelerated by an extrapolation method, adapted to their behavior. 

For periodic chains, the convergence of the gap has been observed to be exponential. In a 
theory with a gap, we expect of course exponential convergence towards the thermodynamic 
limit for a closed chain. This has been shown explicitly in the large-N limit of the nonlinear 



sigma model [15]. We analyze the estimated decay length £(N) at the index N as given by 



P = 2/h p^ , (4) 

\ O-N-2 — a N ) 

where represents the sequence G(N) or e(N). If = A + b exp(— AT/£), then £(N) = £ 
exactly for all N. If a N = A + b/(N - n ) u , then £(N) ~ (N - n )/{v + 1) for N large, and 
the exponent v is given by the asymptotic slope of £(N). 

In Fig. |I|, we have plotted £(iV) for G(N) and e(N). Both curves are concave: the 
estimated exponent v increases with N (for N = 22, respectively, v « 15 and 11). This 
figure shows that the gap and the energy per site converge faster than a power law. This 
is good evidence for the expected exponential behavior of the Haldane conjecture. For this 



reason, we extrapolate with the Shanks transformation p3p^| . We explain in some detail 
this method because we will use it in a different way than Ref. [15| or [25]. If the sequence 
ajy is a sum of k exponentials: 

a N = A + b 1 e~ N ^ + ... + b k e- N ^, (5) 

the limit A is one of the 2k + 1 unknowns and can be obtained by solving the system (J5|) for 
OAr-2fc, ■ ■ ■ , a N, O'N+2, ■ ■ ■ , a N+2k- We call this solution, i.e. the limit A extracted from 
a N-2k, ■ ■ ■ , a N+2k supposing that the sequence is a sum of k exponentials. Of course, if 
the sequence has not exactly this form, then A$ varies with N. The simplest case of 
the Shanks transformation is k — 1. The solution is 

2 

^(1) _ QjV+2 QjV-2 ~ «jy ^ 
a N+2 ~ 2(3jv + O-N-2 

which is also called Aitken's A 2 process. This Aitken-Shanks transformation (^) with k = 1 

*-N 5 



can be iterated [|TJ,^4],^]: it is first applied to a N} then again to the obtained sequence Ay 
and so forth. The iteration of the (biotransformation does not give the A$ , for k > 1. 
For k > 1, the A$ can be computed by using the recursive "cross rule" due to Wynn p4]j26| 



(4v fe+1) - + «(^r 1} - A^)- 1 = (<l 2 - A^)- 1 + (A%> +2 - A%>)-\ (7) 

with the initial conditions A$ = a^, A^ N ^ = oo and where a = 1 (for the Shanks trans- 
formation). The table of the A$ verify many algebraic properties ||24||; for example, if the 
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a,N are the partial sums of a power series, the A$ are the Pade approximants of the series. 
Due to the non-linearity of Eq. (|7|), only a few exact results are known. For example, if 
is a sum of exponentials 



aw 



A + ^Tbi e~ N/il when iV -> oo (8) 



i=l 



with 6 > 6 > " ' " > 0, then f24[ for k fixed and A" — > oo 



A^ ~A + b k+1 (Afc fc +1 ~ Al)2 - • ^±1 ~ A f e -^ +1 with A, = e-^, (9) 

Each /c-iteration removes a exponential and each column k will be more rapidly convergent 
than the previous ones when N goes to the infinity. Even if the exact results are rare, in 
practice a quite general class of sequence is accelerated. 

The parameter a in Eq. (|7|) was introduced by Vanden Broeck and Schwartz [E7| (VBS). 
The table changes when a varies. The Pade-Shanks transformation (at order k) is given 
by a = 1. The iterated Aitken-Shanks (biotransformation is given by a = 0. When 
a = — 1, Hamer and Barber |28j have shown that, if has exactly the power law behavior 
a N = A + 6(1 - i//2)(l - i//4) ... (1 - I//JV) (for AT large, ~ A + V /N v ), the second 

(2) 

column = A for all AT. Hamer and Barber's transformation can be iterated with 
otk = 0, —1, 0, — 1, . . . successively for each column k. 



IV. EXTRAPOLATED VALUES 

On Table O, we give the results of the Shanks transformation (a = 1 and k — 1 to 
5) for the gap G^AQ = E\ — E . The estimated decay lengths £(N, k), defined by Eq. (|J) 
are calculated, for each k, with the A^ 1 . One remarks that the data of each column are 
monotonic. This is also true with the last oblique row. In particular, the £(N,k) decrease 
with k. In fact because of the rather small values of N, £(N, k) represent only an effective 
decay length and £(N, k) involve data from smaller value of A^ than £(N, k — 1). Thus what 
we have to expect is 

C( N , k) < £(N + 2,k- 1). (10) 

These inequalities are all satisfied and the idea that the Shanks transformation removes at 
each step an exponential transient is coherent. The values for N = 2 and 4 do not disturb 
the table. One must compute this table with at least double precision arithmetic (64 bits), 
because large cancellations occur in the cross rule (0). We verified that the precision is 
better than 10 -7 (only 6 digits are retained). 

If the parameter a of Eq. (^) takes a value sufficiently different from a = 1, this table 
loses its consistency (for example, Eq. (|ToD is not verified for all (AT, k) values) and the ex- 
trapolation cannot be trusted. In particular, a = (iterated Aitken-Shanks transformation) 
and at = 0, —1,0, —1, . . . (Hamer and Barber's transformation) do not give a satisfactory 
table. For the gap, the acceptable interval is 0.7 < a < 1.05, (close to a = 1, the pure Shanks 
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transformation) and the respective extrapolations for these bounds are 0.410478 and 
0.410504. Then, with the VBS or Shanks transformation, for periodic chains with N up to 
22, the gap value is 

G(oo) = 0.41049(2). (11) 

The error bar (2.10 -5 ) is an estimation of the systematic error due to the arbitrary parameter 
a. The Shanks transformation does not give a direct estimation of the systematic error with 
respect to the true limit. But the main argument in favor of this method is the good 
regularity of the data of the Table f|. 

Now we study the ground state energy per site e(N) = Eq/N with a same method. In 



Table pi] , we give the results of the Shanks transformation (a = 1 and k = 1 to 5). The 
convergence is faster (in particular, the decay lengths £ are shorter) and the precision is 
better (10 -9 but only 6 digits are given here). For a = 1 (the pure Shanks transformation), 
the table is consistent: all the columns are monotonic, as is the last oblique row, and the 
£'s decrease with k. For a, the acceptable interval is 0.5 < a < 1.1. The respective 
extrapolation for these bounds are —1.401484 and —1.401487. Then, with the VBS or 
Shanks transformation, for periodic chains with N up to 22, the ground state energy per 
site is 

e(oo) = -1.401485(2). (12) 



V. COMPARISONS 

We compare our extrapolations with those published by other authors. For the most 
part, numerical results have been obtained by exact diagonalizations (with simple itera- 
tions or Lanczos method), Monte-Carlo simulations or the real-space renormalization-group 
method. Of course, for a fixed chain length N, we find the same results as the other exact 
diagonalization studies: Blote |4[] for N = 10, Glaus and Schneider || and Parkinson and 
Bonner [7j for iV = 14, Moreo M for N = 16 and Lin M for N = 18. Sakai and Takahashi 



25], with N < 16 results extrapolated by the Aitken-Shanks iterated transformation (Eq. 
D give G(oo) = 0.411(1), compatible with Eq. (0). 

The results of Monte-Carlo calculations have statistical as well as systematic errors but 
longer chains (for example N = 32 or 64) can be studied. The stochastic iteration of Nightin- 
gale and Blote pf gives compatible results: G(32) = 0.413(7) and e(32) = -1.40155(16), as 



the method of Liang [14]: e(64) = —1.402(1). On the other hand, the projector Monte-Carlo 



method of Takahashi JXXJ] is not compatible: e(32) = —1.4023(1). The methods based of 
the Trotter- Suzuki decomposition are characterized by imprecisions when the temperature 
goes to zero and the properties of the ground state are not well reproduced. For example 



Nomura |13| gives G = 0.425. 

Concerning the real-space renormalization-group method (or truncated basis expansion), 
we are in disagreement with some of the published values: Pan and Chen |29j (G = 0.368166 



and e = -1.449724), Mattis and Pan |0J (e = -1.388), and Xiang and Gehring |3J 



(e(oo) = —1.377). On the other hand, we aggre with Lin and Pan JET!: e(oo) = —1.4021(5) 
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and G(oo) = 0.4097(5) and with the recent method of White ]T8| which gives e(oo) = 
-1.401484038971(4) and G(oo) = 0.41050(2). These values are in good agreement with 
ours; the ground state energy is more precise and the precision on the gap value is similar. 
The fact that both methods give results with 6 and 5 identical digits is a good argument 
that there are both quite accurate. 

VI. CORRELATION FUNCTIONS 

In this section, we present our results for the correlation functions 

C N (r) = (-l) r (S* S z r ) = (-iy(SS)/3 , (13) 

for the ground state of a isotropic and periodic chain of length N for iV < 22. Numerical 
values are given in Table [IV]. To compute these quantities, the eigenvalue is not sufficient 
and the eigenvector is required. So the precision for the Cjv(r) is less than for the energies. 
It can be estimated around 10~ 8 , for example by comparing Cjv(1) with E(N), or by direct 
applications of the matrix H on the Lanczos result. 

For an infinite chain, the Haldane conjecture predicts an exponential decrease C^r) ~ 
6exp(— r/£)/y/r when r is large. This is because the underlying continuum theory is a 
nonlinear sigma model which is relativistic in 1+1 dimensions. In fact, if we approximate 
the nonlinear sigma model by free massive bosons, the propagator is the modified Bessel 
function Kq which has this asymptotic behavior. The Haldane conjecture does not deal with 
short-distance details so that there is no a priori preferred choice when trying to fit data on 
the full range of spin-spin separation. 

For a periodic chain, one has the equality CV(r) = C^(N — r). To extract the correlation 
length £ for N and r large, we analyze our data with the guess C^v(r) = C^ir) + C OQ (N — r), 
which is reasonable if £ -C N. For some classical spin systems (one-dimensional Potts model, 
Ising chain with a magnetic field, . . . ) with non-vanishing temperature, the corrections to 
this formula are of order 0(Coo(N)). 

We have verified that exp(— r/£)/ \fr fits the data better than exp(— r/£) or exp(— r/£)/r. 
From these three forms, the estimated values for £ are respectively 6.2, 4.5 and 10, for 
N = 22. In Fig. |2], we compare the Bessel function Ko(r/£) and exp(— r/£)/y/r. Both fits 
are comparable in quality. But the estimated correlation lengths £ are slightly different: 5.9 
versus 6.2 for N = 22. For iV = 10, they are respectively 5.4 and 6.2. We notice that the 
optimal £ is 6.2, for all iV < 22, with exp(— r/£)/-y/r, whereas the estimation for £ with the 
Bessel function vary with N. For this reason, we prefer the former but we keep in mind that 
both laws have the same asymptotic behavior and that only a exact solution of the model 
can give the correlations for short range. 

It is interesting to compare the correlation length £ (obtained with Cjv(r)) and the decay 
length £(iV) (Eq. |) of the gap G(N) (Table 0) and energy e(N) (Table ^j. On Fig. 
[TL we see that it is not excluded that £ = 6.2 is the limit of £(iV). The extrapolation of 
£(iV) with the Shanks transformation gives 5.5 for the gap and 4.6 for the energy, but the 
columns of the Shanks table are non-monotonic and these results are only qualitative. Of 
course, by analyzing the convergence of the G(N) and e(N) with an exponential corrected 
by a power law (as 1/vN), the estimations of £(iV) are greater, and the extrapolations are 
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closer from 6.2. Then this comparison is only qualitative and requires longer chains. Our 
value (£ = 6.2) is equal to the estimate quoted by Nomura [13| (£ = 6.25) and by Liang 



jnj (£ = 6.2) with Monte Carlo methods for A = 64. It is comparable with the results of 
Takahashi pf (f = 5.5 ±2.) with a Monte Carlo method for A = 64, by Kubo || (£ = 6.7) 
with a transfer-matrix method and by White and Huse |18j (£ = 6.03) with the real-space 
renormalization-group method. 



VII. CONCLUSION 

The main limitation of exact diagonalizations is, of course, the small lengths that can 
be studied. The numerical complexity grows as the exponential (S(S + 1))^ for iV spins S 
and the limits of computer power are fastly reached. The length A of the system must be 
compared with the physical correlation length £, and in fact, the situation for the S=l AF 
spin chain is quite favorable. Within the Haldane conjecture, £ is finite for integer spins and 
shortest for small S. We have shown that some quantities can be measured with excellent 
accuracy: the gap and the ground state energy. On the other hand, the correlations Coo(r) 
(and thus the correlation length £) clearly need longer chains. 

The main advantages of exact diagonalizations are that they depend only on one pa- 
rameter (the size N) and give exact results (i.e. with machine precision). One has to deal 
only with the thermodynamic limit. By comparison, methods based on Trotter-Suzuki de- 
composition have three parameters (number of slices, temperature and length of the chain) 
and systematic errors which decrease by extrapolating in the number of slices. Monte-Carlo 
methods have their own parameters (number of walkers or length of simulations, . . . ) which 
must be tuned, and the results have statistical fluctuations as well as systematic errors. 
Real-space renormalization-group methods have to extrapolate w.r.t. the number of basis 
states and the chain length. 

The high precision allows the use of sophisticated extrapolation methods and we are able 
to validate some assumptions on the asymptotic behavior. Fig. [I] suggests that the use of 
the Shanks transformation is optimal concerning gap extrapolation. In fact the parameter 
a of the more general VBS transformation can vary only in a small interval around a — 1. 
This shows that our choice is not arbitrary but dictated by the data. The results of exact 
diagonalizations combined with a careful extrapolation can give physical quantities in the 
thermodynamic limit with a good precision. 



APPENDIX: PROGRAMMING TECHNIQUES 

We used a Cray 2 of the CEA with a central memory of 256 Megawords of 64 bits. Some 
details of our program are useful only for this kind of machine in particular and are not 
described here. The algorithm has two main parts, the building of the sparse matrix H (and 
its storage on disks) and the matrix-vector H ■ V multiplication, needed for the Lanczos 
iterations. 

We consider first the matrix multiplication. The matrix H is very sparse. For A = 22, 
its order is ~ 37 million and the number of non- vanishing elements is (on average) 8 A/9 per 
rows (when A is large). We use the classical storage by rows with only the non- vanishing 
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elements of H (values and column number) stored. In practice, one bit is needed for the 
value (±1) and 26 bits for the column number. So, two elements are stored in a 64-bits word 
and 3.5 Gigabytes are used for H . The matrix-vector multiplication is done by an indirect 
addressing of the elements of the vector, where the address is the column number. This 
indirect addressing is the most time consuming part of the program and it is intrinsic to this 
sparse storage method. For N = 22, a multiplication needs 190 seconds of one Cray 2 CPU. 

The most difficult part of the algorithm is the building of the matrix with use of symme- 
tries. Each z-axis basis state \s\, . . . , s^) is described by the number J2r( s r + l)3 r_1 . First, 
the list of the symmetrized basis states is obtained. Each symmetrized state is represented 
by the state of the z-axis basis, which contributes and has the smallest number. Then the 
hamiltonian H operates on this list and generates other states. The problem is to find the 
corresponding symmetrized states (and their phases). Possible methods are a) each gener- 
ated state is symmetrized by action of all the symmetries operators or b) a storage table 
gives, for each 2-axis state, the symmetrized one and the generated state is searched in this 
table. The first method is too time consuming and the second one uses too much memory. 

We use an intermediate method with a decomposition in two sublattices ||, A = {s2 r } 
and B = {s2r+i}- The symmetries R x , Lr and T 2k do not exchange A and B. We call them 
sublattice symmetries. On the other hand, symmetries T 2k+l exchange A and B. Then, 
each symmetry is the product of a sublattice symmetry and possibly T. Since a sublattice 
is described only by 3 N / 2 states, we can use a storage table which gives for each sublattice 
state the symmetrized one. For all the chain, A-symmetrization consist to symmetrized A 
and to operate on B with the same operator. Since a storage table of size 3 N ^ 2 is available, 
it does not require much time or memory. The last step is the action of T, which exchange 
A and B, and S z , which imposes S Z (A) + S Z (B) = 0. By symmetrizing by T, the number of 
the A-symmetrized states (around 78 millions for N = 22 and S z = 0) is divided by two (for 
N large). On our computer, we keep on memory the list which gives the fully symmetrized 
state for each A-symmetrized one. This list has some properties of factorization, as well 
explained in Ref. 0, and the location of each state is easily obtained by considering each 
sublattice. To summarize, a generated state is, in a first step, symmetrized by R x , Lr and 
T 2k (which let invariant the sublattices), and in a last step by T. The first step needs only 
short lists (3 N ^ 2 ) and the final one a big list, for which the length is two times the order of 
the fully symmetrized block. For iV = 22, our program needs 2200 seconds of one Cray 2 
CPU to built the matrix H. 

For one block (|0) or |1, 0)), with the Lanczos method, the precision cannot be improved 
after 55 iterations (for N = 22). To compute eigenvectors, the Lanczos method is not 
optimal because all the intermediate vectors must be stored. To do that, 16 Gigabytes are 
required. Then, a first Lanczos calculation gives the eigenvalues and the coordinates of 
eigenvectors on the Lanczos basis. A second Lanczos calculation is needed to generate the 
eigenvectors. 

The computations of this article have used 12 hours of one Cray 2 CPU. 
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FIGURES 



FIG. 1. The estimated decay length £(iV) given by Eq. (Q) for the gap G(N) and the ground 
state energy per site e(N) vs. the length N of the periodic chain. 

FIG. 2. The ratio of the correlation function CV(r) = (— l) r (Sq ■ S*) and two proposed laws 
versus r. The CV(r) have been exactly computed for N = 22. The ratios are normalized to 1 for 
r = N/2. 
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TABLES 



TABLE I. Dimension of the largest block (S z = 0,k = 0, Lr = +1,R X = +1), ground state 
energy Eq, first excitation energy E\, gap G(N) = E\ — Eq and ground state energy per site 
e(N) = Eq/N for chain lengths A = 2 to 22. These results are obtained by exact diagonalization. 
Previous results for N = 8 are given by De Neef @, A = 10 by Blote [{§, N = 12 by Botet and 
Jullien ||, = 14 by Glaus and Schneider || and Parkinson and Bonner 0, N = 16 by Moreo 
|g and N = 18 by Lin @. 



N 


dimension 


—Eq 


-E x 


gap G(N) 


-e(AT) 


2 




2 


4.0 


2.0 


2.0 


2.0 


4 




5 


6.0 


5.0 


1.0 


1.5 


6 




15 


8.617423181814 


7.896795819190 


0.720627362624 


1.436237196969 


8 




59 


11.336956077897 


10.743400823522 


0.593555254375 


1.417119509737 


10 




290 


14.094129954932 


13.569322004518 


0.524807950414 


1.409412995493 


12 


1 


728 


16.869556139477 


16.385359669563 


0.484196469914 


1.405796344956 


14 


11 


549 


19.655133499935 


19.196168152997 


0.458965346938 


1.403938107138 


16 


82 


790 


22.446807281171 


22.004011719811 


0.442795561360 


1.402925455073 


18 


617 


898 


25.242312007671 


24.810090537803 


0.432221469868 


1.402350667093 


20 


4 730 


966 


28.040291720480 


27.615081406019 


0.425210314461 


1.402014586024 


22 


36 871 


567 


30.839898879910 


30.419383859516 


0.420515020394 


1.401813585450 



TABLE II. The Shanks extrapolations Afi' for the gap values G(N) = E 1 - E with k = 1 to 
5 (with the VBS parameter a = 1). The estimated decay lengths £ are obtained by applying the 



formula @ on the A^ k x) . 



N 


G(N) 








e 


A® 




AW 






2 


2.000000 




















4 


1.000000 


0.612320 


1.57 
















6 


0.720627 


0.487533 


2.54 


0.435259 


1.91 












8 


0.593555 


0.443776 


3.26 


0.417985 


2.28 


0.412584 


1.59 








10 


0.524808 


0.425578 


3.80 


0.413089 


2.43 


0.411146 


1.75 


0.410712 


1.57 




12 


0.484196 


0.417574 


4.20 


0.411524 


2.53 


0.410744 


1.98 


0.410555 


2.08 


0.410498 1.85 


14 


0.458965 


0.413941 


4.50 


0.410954 


2.64 


0.410590 


2.32 


0.410501 


2.38 




16 


0.442796 


0.412240 


4.71 


0.410714 


2.77 


0.410523 


2.69 








18 


0.432221 


0.411414 


4.87 


0.410600 


2.94 












20 


0.425210 


0.410996 


4.99 
















22 


0.420515 
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TABLE III. The Shanks extrapolations A N for the ground state energy per site e(N) = Eq/N 
with k = 1 to 5 (with the VBS parameter a = 1). The estimated decay lengths £ are obtained by 
applying the formula (ffl) on the A^ -1 ). 



AT 


-e(N) 


A w 


£ 


AW 


e 


A(3) 


£ 


AW 


£ 


A(5) 


£ 


2 


2.000000 






















4 


1.500000 


1.426917 


0.97 


















6 


1.436237 


1.408933 


1.66 


1.403743 


1.50 














8 


1.417120 


1.404208 


2.20 


1.402154 


1.86 


1.401683 


1.56 










10 


1.409413 


1.402598 


2.64 


1.401712 


2.11 


1.401544 


1.75 


1.401503 


1.59 






12 


1.405796 


1.401974 


3.00 


1.401571 


2.30 


1.401505 


1.95 


1.401490 


1.83 


1.401486 


1.68 


14 


1.403938 


1.401713 


3.29 


1.401521 


2.47 


1.401492 


2.17 


1.401486 


2.01 






16 


1.402925 


1.401596 


3.53 


1.401500 


2.64 


1.401487 


2.38 










18 


1.402351 


1.401541 


3.73 


1.401492 


2.82 














20 


1.402015 


1.401514 


3.89 


















22 


1.401814 























TABLE IV. The correlation functions Cjv(r), calculated by exact diagonalization for up to 
22. For < 18, these results have been published by other authors @-9 



r 


N = 6 


8 


10 


12 


14 


16 


18 


20 


22 


1 


0.47874573 


0.47237317 


0.46980432 


0.46859878 


0.46797936 


0.46764181 


0.46745022 


0.46733819 


0.46727118 


2 


0.28844542 


0.27210249 


0.26392567 


0.25918542 


0.25622175 


0.25429251 


0.25300867 


0.25214355 


0.25155626 


3 


0.28606604 


0.24086913 


0.22135314 


0.21075706 


0.20436261 


0.20028789 


0.19761348 


0.19582835 


0.19462472 


1 




0.21561295 


0.18479849 


0.16814782 


0.15810415 


0.15169940 


0.14749166 


0.14468055 


0.14278366 


5 






0.18180007 


0.15402811 


0.13824755 


0.12849332 


0.12220081 


0.11804519 


0.11526281 


6 








0.14543474 


0.12353893 


0.11017252 


0.10161828 


0.09599964 


0.09225227 


7 










0.12121726 


0.10228686 


0.09052395 


0.08293057 


0.07792178 


8 












0.09842421 


0.08305179 


0.07322701 


0.06679153 


9 














0.08143053 


0.06847270 


0.06012605 


10 
















0.06646187 


0.05588756 


11 


















0.05479614 
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